%%%-------------------------------------------------------------------
%%% File    : p46.erl
%%% Author  :  Plamen Dragozov <plamen at dragozov.com>
%%% Description : 
%%% It was proposed by Christian Goldbach that every odd composite 
%%% number can be written as the sum of a prime and twice a square.
%%%
%%% 9 = 7 + 2*1^(2)
%%% 15 = 7 + 2*2^(2)
%%% 21 = 3 + 2*3^(2)
%%% 25 = 7 + 2*3^(2)
%%% 27 = 19 + 2*2^(2)
%%% 33 = 31 + 2*1^(2)
%%%
%%% It turns out that the conjecture was false.
%%%
%%% What is the smallest odd composite that cannot be written as the sum of a prime and twice a square?
%%%
%%%
%%% Created :  1 Jan 2009
%%%-------------------------------------------------------------------
-module(p46).

%% API
-compile(export_all).

% X = 2*n + 1
% X not prime
% X = p + 2*a**2, p - prime, a - integer

%%====================================================================
%% API
%%====================================================================
%%--------------------------------------------------------------------
%% Function: 
%% Description:
%%--------------------------------------------------------------------
solution() ->
    for(9, [2, 3, 5, 7]).

%%====================================================================
%% Internal functions
%%====================================================================
% For every odd X
% -> If it is not prime:
% iterate through the primes < X and check if it meets the condition 
% for any of them
% -> If it is prime add to the list of primes
for(X, Primes) ->
    case is_prime(X, math:sqrt(X), Primes) of
        true ->
            for(X + 2, Primes ++ [X]);
        _ -> case test(X, Primes) of
                 true ->
                     for(X + 2, Primes);
                 _ ->
                     X
             end
    end.

test(_, []) -> false;
test(X, [H|T]) ->
    A = math:sqrt(((X - H) / 2)),
    case A > trunc(A) of
        true ->
            test(X, T);
        _ -> true
    end.

is_prime(_, _, []) -> true;
is_prime(X, Sqrt, [H|T]) ->
    X rem  H =/= 0 andalso (H > Sqrt orelse is_prime(X, Sqrt, T)).

